%Figure 3 - ROC Curves%
%Franzese, Hays, Cook 2014%

clear;

load 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM Conditional Accept\Replication Folder\Application - Civil War\Results\hb_ksg_RR_rho_phi.mat'

pt=pt';

mult_fx = inv(i_nt - pt(1,10)*TL - pt(1,11)*WNT);
SIGMA = mult_fx*mult_fx';
rcf = mult_fx*(pt(1,1) + pt(1,2)*XNT(:,1) + pt(1,3)*XNT(:,2) + pt(1,4)*XNT(:,3)+ pt(1,5)*XNT(:,4)+ pt(1,6)*XNT(:,5)+ pt(1,7)*XNT(:,6)+ pt(1,8)*XNT(:,7)+ pt(1,9)*XNT(:,8));

nt = 1434;
n = 31;

y = Y;

pred_prob = zeros((nt-n),1); 

for i = n+1:nt; 
    ind_TL = find(TL(i,:)); 
    ind_TL_2 = find(TL(ind_TL,:)); 
    ind_TL_3 = find(TL(ind_TL_2,:));
    %ind_TL_4 = find(TL(ind_TL_3,:));

    ind_W = find(WNT(i,:));
    %id = [ind_TL ind_TL_2 ind_W]; 
    id = [ind_TL ind_TL_2 ind_TL_3 ind_W];
    %id = [ind_TL ind_TL_2 ind_TL_3 ind_TL_4 ind_W]; 
    id_outcome = y(id); 
    
    id_1 = id(id_outcome~=0);
    id_0 = id(id_outcome==0);
    
    ub = Inf(1,numel(id_0)); 
    lb = -Inf(1,numel(id_1));
        
    p_B = mvncdf([lb, rcf(id_0)'], [rcf(id_1)', ub],0,SIGMA([id_1,id_0],[id_1,id_0])); 

    p_AB = mvncdf([lb, rcf(id_0)',-Inf], [rcf(id_1)', ub, rcf(i)],0,SIGMA([id_1,id_0,i],[id_1,id_0,i])); 

    pred_prob(i-n) = p_AB/p_B;         
end 

%Plot ROC curves 
[X,Y,T,AUC_sprob] = perfcurve(y(n+1:nt,1), pred_prob, 1);

ytl = TL*y(:,1); 
ysl = WNT*y(:,1); 

y = y(n+1:nt,:); 
XNT1 = XNT(n+1:nt,:);
ytl1 = ytl(n+1:nt,:);
ysl1 = ysl(n+1:nt,:); 
py1 = py(n+1:nt,:);
cons = ones((nt-n),1); 

X_mat = [cons XNT1 ytl1 ysl1]; 
b_rs = glmfit(X_mat,y,'binomial','link','probit','constant','off');
p_rs = glmval(b_rs,X_mat,'probit','constant','off');

[X,Y,T,AUC_rs] = perfcurve(y, p_rs, 1);


X_dth = [cons XNT1 py1 ysl1]; 
b_dth = glmfit(X_dth,y,'binomial','link','probit','constant','off');
p_dth = glmval(b_dth,X_dth,'probit','constant','off');

[X,Y,T,AUC_dth] = perfcurve(y, p_dth, 1);

AUC_total = [AUC_sprob AUC_rs AUC_dth]

y_conflict = y(ytl1==1); 
y_peace = y(ytl1==0); 

[X_stp_peace,Y_stp_peace,T,AUC_stp_peace] = perfcurve(y_peace, pred_prob(ytl1==0), 1);
[X_rs_peace,Y_rs_peace,T,AUC_rs_peace] = perfcurve(y_peace, p_rs(ytl1==0), 1);
[X_dth_peace,Y_dth_peace,T,AUC_dth_peace] = perfcurve(y_peace, p_dth(ytl1==0), 1);

AUC_peace = [AUC_stp_peace AUC_rs_peace AUC_dth_peace]

subplot(1,2,1)
plot(X_stp_peace,Y_stp_peace,'-',X_rs_peace,Y_rs_peace,'--',X_dth_peace,Y_dth_peace,':')
xlabel('False positive rate'); ylabel('True positive rate')
title('ROC when observed peace at t-1')

[X_stp_con,Y_stp_con,T,AUC_stp_con] = perfcurve(y_conflict, pred_prob(ytl1==1), 1);
[X_rs_con,Y_rs_con,T,AUC_rs_con] = perfcurve(y_conflict, p_rs(ytl1==1), 1);
[X_dth_con,Y_dth_con,T,AUC_dth_con] = perfcurve(y_conflict, p_dth(ytl1==1), 1);

AUC_conflict = [AUC_stp_con AUC_rs_con AUC_dth_con]

subplot(1,2,2)
plot(X_stp_con,Y_stp_con,'-',X_rs_con,Y_rs_con,'--',X_dth_con,Y_dth_con,':')
xlabel('False positive rate'); ylabel('True positive rate')
title('ROC when observed conflict at t-1')
legend('STP-MSL','RS-ML','DTH-ML','Location','southeast')


